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Abstract 

An asynchronous, variational method for simulating elastica in complex contact and impact scenarios 
is developed. Asynchronous Variational Integrators [Tj (AVIs) are extended to handle contact forces 
by associating different time steps to forces instead of to spatial elements. By discretizing a barrier 
potential by an infinite sum of nested quadratic potentials, these extended AVIs are used to resolve con- 
qv tact while obeying momentum- and energy-conservation laws. A series of two- and three-dimensional 

i—i examples illustrate the robustness and good energy behavior of the method. 
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1. Introduction 
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Variational integrators (Vis) 2, 3, 4] are a general class of time integration methods for Hamiltonian 
systems whose construction guarantees certain properties highly desirable of numerical simulations. 
Instead of directly discretizing the smooth equations of motion of a system, the variational approach 
instead discretizes the system's action integral. By analogy to Hamilton's least action principle, a 
discrete action can be formed, and discrete Euler-Lagrange equations derived by examining paths 
which extremize it. From the Euler-Lagrange equations, discrete equations of motion are readily 
recovered. As a consequence of this special construction, Vis are guaranteed to satisfy a discrete 
formulation of Noether's Theorem [5], and as a special case conserve linear and angular momentum. 
Vis are automatically symplectic [6] ; while they do not necessarily conserve energy, conservation of the 
symplectic form assures no-drift conservation of energy over exponentially many time steps [6] . 

Given the many advantages of Vis, it is natural to apply them to the handling of contact and impact, 
a long-studied and challenging problem in physical simulation. Unfortunately, a naive application of a 
• • contact algorithm to a variational integrator is not guaranteed to preserve the variational structure of 

the time integration method, and in practice one observes that the good energy behavior is lost. For 
this reason, a few recent works have explored structure-preserving approaches for contact mechanics [7j 
[SI US]- Common to all these approaches is a synchronous treatment of global time, in which the 
entire configuration is advanced from one intant in time to the next. While synchronous integration 
is attractive for its simplicity, it has the drawback that a spatially-localized stiff mode — such as that 
associated to a localized contact — can force the global configuration to advance at fine time steps. 

Indeed, mechanical systems are almost never uniformly stiff. Different potentials have different 
stable time step requirements, and even for identical potentials this requirement depends on element 
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size, since finer elements can support higher-energy modes than coarser elements. Any global time 
integration scheme cannot take advantage of this variability, and instead must integrate the entire 
system at the globally stiffest time step. Suppose the system can be partitioned into elements such 
that each force acts entirely within one element. Then asynchronous variational integrators (AVIs) pQ 
generalize Vis by allowing each element to have its own, independent time step. Coarser elements can 
then be assigned a slower "clock," and finer elements a faster one. Asynchrony avoids the undesirable 
situation in which a small number of very fine elements degrade overall performance. AVIs retain 
all of the properties of variational integrators mentioned above, except for discrete symplecticity. 
However, AVIs instead preserve an analogous discrete multisymplectic form, and it has been shown 
experimentally that preservation of this form likely induces the same long-time good energy behavior 
that characterize symplectic integrators pQ. 

To our knowledge, this work is the first to consider an asynchronous, variational treatment of contact 
shown to retain multisymplecity. We are also aware that Ryckman and Lew |11| are concurrently 
investigating extending the AVI framework to incorporate contact response. 

The starting point for this approach is the selection of the penalty method as a model for contact [121 
113) . For each pair of elements in the system, a potential is added that is (piecewise) quadratic in the gap 
function measuring the separation distance between the two elements. This potential vanishes when 
elements are sufficiently far apart, and increases with increasing interpenetration, so that approaching 
elements feel a force that resists impact. This approach suffers two limitations, however. Firstly, these 
contact potentials are fundamentally nonlocal phenomena: for every pair of elements that might come 
into contact during the course of the simulation, a potential coupling the two must be added. As will 
be shown, the fact that contact potentials cannot be expressed as the integration over the material 
domain of an energy density depending only on a neighborhood of the domain will present a technical 
obstruction to the original formulation of AVIs, but fortunately one that can be overcome by a natural 
generalization. 

Secondly, penalty forces have a well-studied performance-robustness tradeoff [13]: adding a half- 
quadratic potential requires choosing an arbitrary stiffness parameter, and for any stiffness chosen 
for the penalty potential, two approaching elements will interpenetrate some distance, and in the 
worst case tunnel completely through each other. Moreover, the stable time step of the penalty force 
decreases as stiffness increases, so choosing a very stiff penalty potential is untenable as a solution 
to excessive penetration or tunneling. In practice, users of the method must determine an adequate 
penalty stiffness by iterated tweaking of parameters, until the simulation completes without collision 
artifacts. An appealing modification of the penalty approach replaces the quadratic potential with 
a nonlinear barrier potential [T5] that diverges as the configuration approaches contact. Because the 
barrier diverges, its stiffness is unbounded, necessitating a time-adaptive time stepping method. This 
work presents a discrete analogue of the barrier potential — an infinite sequence of discrete penalty 
layers — that in effect enables AVIs to serve as adaptive integrators. 

This paper 

• extends the construction of AVIs so that a discretization into disjoint elements is no longer 
necessary, by associating a clock to each force instead of to each element; 

• demonstrates that this generalization does not destroy the desirable integration properties guar- 
anteed by the variational paradigm, most importantly the conservation of the discrete multisym- 
plectic form; 

• leverages this extension to equip the AVI framework with a contact model. The proposed barrier 



method uses a divergent sequence of quadratic potentials that guarantees non-penetration and 
retains the asynchrony or conservation properties of AVIs; 

• presents numerical evidence to support the claim that by retaining the symplectic structure 
of the smooth system, simulations of thin shells undergoing complex (self-)interactions have 
demonstrably good long-time energy and momentum behavior; 

• describes simple extensions to the contact model to allow for controlled, dissipative phenom- 
ena, such as a coefficient of restitution and kinetic friction. Although there is not yet theory 
explaining the energy behavior of dissipative simulations run under a variational integrator, em- 
pirical evidence is presented to show that the proposed method produces smooth, controlled, and 
qualitatively correct energy decay. 

This paper complements the publication |16j , which provides a detailed description of the software 
implementation using kinetic data structures |17j . For completeness, Section p] briefly introduces these 
concepts. 

2. Related Work 

The simplest contact models for finite element simulation follow the early analytical work of 
Hertz |18j in assuming frictionless contact of planar (or nearly planar) surfaces with small strain. 
In this regime, several approaches have been explored to arrive at a weak formulation of contact; 
for a high-level survey of these approaches, see for example the overview by Belytschko et al. [19] or 
Wriggers 20J. The first of these are the use of penalty forces, described for instance by Oden [5T] and 
Kikuchi and Oden |22j . The penalty approach results in a contact force proportional to an arbitrary 
penalty stiffness parameter and to the rate of interpenetration, or in more general formulations to an 
arbitrary function of rate of interpenetration and interpenetration depth; Belytschko and Neal 14J 
discuss the choosing of this parameter in Section 8. Recent work by Belytschko et al. [23] uses moving 
least squares to construct an implicit smooth contact surface, from which the interpenetration distance 
is evaluated. Peric and Owen |24| describe how to equip penalty forces with a Coulomb friction model. 

Seeking to exactly enforce non-penetration along the contact surface leads to generalizations of 
the method of Lagrange multipliers. Hughes et al. [35] and Nour-Omid and Wriggers [5B] provide 
an overview of this approach in the context of contact response. Such contraint enforcement can 
be viewed as a penalty force in the limit of infinite stiffness, impossible to attain in practice since 
the system becomes ill-conditioned. Taylor and Papadopoulos [27J considers persistent contact by 
extending Newmark to treat jump conditions in kinematic fields, thus reducing undesirable oscillatory 
modes. However, the effects of these modifications on numerical dissipation and long-time energy 
behavior is not considered. 

The Augmented Lagragian method blends the penalty and Lagrange multiplier approaches, and 
combines the advantages of both: unlike for pure penalty forces, convergence to the exact interpene- 
tration constraint does not require taking the penalty stiffness to infinity, and the Lagrange multiplier 
solve tends to be well-conditioned. Bertsekas [55] gives a mathematical overview of the augmented 
Lagrangian method, and Wriggers et al. [29] and Simo and Laursen [30] expand on its application to 
contact problems in finite elements. 

Non-smooth contact requires special consideration, since in the non-smooth regime there is no 
straightforward way of defining a contact normal or penetration distance. Simo et al. |31] discretize 
the contact surface into segments over which they assume constant contact pressure; this formulation 



allows them to handle non-node-to-node contact using a perturbed Lagrangian. Kane et al. |32j apply 
non-smooth analysis to resolve contact constraints between sharp objects. Pandolh ct al. [TU] extend the 
work of Kane et al. by describing a variational model for non-smooth contact with friction. Cirak and 
West |33| decompose contact resolution into an impenetrability-enforcement and momentum-transfer 
step, thereby exactly enforcing non-interpenetration while nearly conserving momentum and energy. 

Several authors have explored a structure-preserving approach to solving the contact problem. 
Barth et al. [7] consider an adaptive-step-size algorithm that preserves the time-reversible symmetry of 
the RATTLE algorithm, and demonstrate an application to an elastic rod interacting with a Lennard- 
Jones potential. Kane et al. [8] show that the Newmark method, for all parameters, is variational, 
and construct two two-step dissipative integrators that yield good energy decay. Laursen and Love [5] , 
by taking into account velocity discontinuities that occur at contact interfaces, develop a momentum- 
and energy-preserving method for simulating frictionless contact. This paper shares with these last 
approaches the viewpoint that structured integration, with its associated conservation guarantees, is 
an invaluable tool for accurately simulating dynamic systems with contact. 

Although several previous approaches are also adaptive, the algorithm described in this paper is 
the first structured integrator for contact mechanics that achieves time adaptivity using asynchrony. 
This novel approach guarantees the robustness of the proposed integrator, without compromising the 
good properties of structured integration. 

3. Variational Integrators 

This section presents a background on variational integration and symplectic structure [6j SI [5] . 

Let j(t) be a piecewise-regular trajectory through configuration space Q, and j(t) = ^{t) be the 
configurational velocity at time t. For simplicity, assume that the kinetic energy of the system T de- 
pends only on configurational velocity, and that the potential energy V depends only on configurational 
position, so that the Lagrangian L at time t may be written as 

L(q,q)=T(q)-V(q). (1) 

Then given the configuration of the system q at time to and qf at tf, Hamilton's principle |34j 
states that the trajectory of the system j(t) joining 7(^0) — 9o and 7(£/) = qf is a stationary point of 
the action functional 

S( 7 )= / '' I[ 7 (t),7(*)]# 

Jt 

with respect to taking variations $7 of 7 which leave 7 fixed at the endpoints to, tf. In other words, 
7 satisfies 

dS(j) • <5 7 = 0. (2) 

Integrating by parts, and using that 5j vanishes at to and t\, 

(^(7,7) ■ 57+^(7,7) • *t) dt = J (--^-(7) " 7^(7)7j • hdt = 0. 



Since this equality must hold for all variations S~/ that fix 7's endpoints, 
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dt \ dq 



(7) = 0, 



(3) 



the Euler- Lagrange equation of the system. This equation is a second-order ordinary differential 
equation, and so has a unique solution 7 given two initial values 7(^0) & n d 7(^0) ■ 

3.1. Symplecticity 

The flow O s : [7(f), 7(i)] H> py(i + s),7(£ + s)] induced by pi) has many structure-preserving prop- 
erties; in particular it is momentum-preserving, energy-preserving, and symplectic |35j . To derive this 
last property, for the remainder of this section the space of trajectories is restricted to those that satisfy 
the Euler-Lagrange equations. For such trajectories, if the requirement that #7 fix the endpoints of 7 
is relaxed, then the boundary terms of the integration by parts are no longer and 



dT 

dS{i) ■ h = -qT [^ q {q,q)\ ■ h 



(4) 



where n q is projection onto the second factor. 

Since initial conditions (q, q) are in bijection with trajectories satisfying the Euler-Lagrange equa- 
tion, such trajectories 7 can be uniquely parametrized by initial conditions [7(^0) > 7(^0)] • For the 
remainder of this section variations ^7 are also restricted to first variations: those variations in 
whose direction 7 continues to satisfy the Euler-Lagrange equations. These are also parametrized 
by variations of the initial conditions, (5q,Sq). For conciseness of notation, the change of variables 
v{t) = (7(f), 7"(t)) and 8v(t) = \5^{t), Sj(t)] can be used; using this notation the above two facts can be 
rewritten as u(t) = &t-t ^(to) and 8v(t) = Qt-t 0if 5v{to). The action 111), a functional on trajectories 
7, can also be rewritten as a function Si of the intial conditions, 



Si(q,q) 



L[e t (q,q)]dt, 



so that 



dS( 7 ) • 8 1 = dSi Ht )} ■ 8v{t ). 
Substituting all of these expressions into Q, 



dT 



dS z [u(to)] ■ 5v(to) = ( jt o ir q ) [9 t _ to t/(*o)] • h(t) 



dT 



o TT q [9 t _ to j/(t )] dq ■ 5v{t) 



dq ' "" 
dT 

{®t s -t a *0L ~ 0£)i/(t o ) ' <M*o), 



o TT q [Q t - to v{t )} dq ■ Q t -t Jv(t ) 



where $l is the one-form I ?? o n^ j dq. Since dSj is exact, 



d*S i = = &t t -t *d6 L -d9L, 

so since to and tf are arbitrary, Q*d0L — dO^ for arbitrary times s, and O preserves the so-called 
symplectic form dOj^ . 

3.2. Discretization 

Discrete mechanics [3JH HJ |371 [351 HI E] describes a discretization of Hamilton's principle, yielding a 
numerical integrator that shares many of the structure-preserving properties of the continuous flow S . 
Consider a discretization of the trajectory 7 : [to, tf] — > Q by a piecewise linear trajectory interpolating 
n points q = {go, Qi, ■ ■ ■ Qn-i}, with qo = 7(^0) an d q n -i — l(tf), where the discrete velocity ft+1/2 on 
the segment between q^ and qi+\ is 

_ Qi+i — Qt u tf — to 
Qi+1/2 — r 7 n — ■ 

An analogue of fcfy in this discrete setting is needed. To that end, a discrete Lagrangian 

L d (q a ,q b )=T(^^j-V(q b ) (5) 

can be formulated, as well as a discrete action 



Sd(l) = ^2 hL d (qi,q i+ i). (6) 



i=0 



Motivated by |2]), a discrete Hamilton's principle can be imposed: 

dS d (q) ■ 8q = 

for all variations 5q = {Sqo, 8q\ 1 . . . , 8q n -\} that fix q at its endpoints, i.e., with 8qo = 5q n -\ = 0. For 
ease of notation, the kinetic and potential energy terms in ^j can be written to depend on (q a ,q b ), 
two points of phase space consecutive in time, instead of (q, q): 

rr , \ r r fQb-q a \ rr> f s dT fq b -q a 

TdiQa, q b ) = T[ T d (q a , q b ) = 



h ) -a^u.,^, g . y h 

dV 

dq 



dV 
V d {qa,q b ) = V(q b ) Vd(la,qb) = -K-(Qb) 



Then 

dS d (q) -6q=^2h (D^fo, q i+1 ) ■ Sq, t + D 2 L d {q l , q l+1 ) ■ Sq i+1 ) 



i=0 

^ h \ ~^ T d(*>ft+i) '<*& + J, T d(QuQi+l) -5q l+ i - — (q i+1 ) -5q l+ i 

dV 
T d (q n -2,q n -i) ■ hn-l - T d (q 0l qi) ■ Sq - h — {q n -i) ■ 5q n -i 

+ 5Z ( T d(9i-i)&) -T d (qi,q i+ i) - h^ 

ra— 2 / ota 



i=l x 

Since <5gi is unconstrained for 1 < i < n — 2, 

<9T <9T dV 

gvWi+1/2) - -^rfe-1/2) = -h-Q-w), i = l,...,n-2, (7) 

the discrete Euler-Langrange equations of the system. 

Unlike in the continuous settings, the discrete Euler-Lagrange equations do not always have a 
unique solution given initial values qo and q\. Therefore in all that follows it is assumed that T d and 
V d are of a form so that |7| gives a unique qi+i given qi and qi-\ — this assumption always holds, for 
instance, in the typical case where T d is quadratic in q. Then the discrete Euler-Lagrange equations 
give a well-defined discrete flow 

F : (qi-i,qi) t-> (quqi+i), 
which recovers the entire trajectory from initial conditions, in perfect analogy to the continuous setting. 

3.3. Symplecticity of the Discrete Flow 

By analogy to the continuous setting, it is desired that F preserve a symplectic form, just as dO^ 
is preserved by 0. As in the continuous setting, trajectories q are restricted to those that satisfy 
the discrete Euler-Lagrange equations, and variations to first variations (and the condition that these 
variations vanish at the endpoints is lifted), yielding 

dV 
dS d {q) ■ Sq = T d (q n _ 2 , q n -i) • Sq n -i - T d (q ,qi) ■ 5q a - h-—{q n _ 1 ) ■ £g n _i. 

F k denotes the discrete flow F composed with itself k times, or k "steps" of F. Again, all q satisfying 
¥l\ can be parametrized by initial conditions vq = ((70,91), and first variations by 5 vq — (6qo,5qi), so 
that the discrete action can be rewritten as 

n-2 

4=0 



Putting together all of the pieces, 
dS id (v ) ■ 8vq = dS d (q) ■ 5q 



dV 
= T' d {q n -2,q n -i) ■ 5q n -x - T d (q , q{) ■ 5q - h—(q n ^ 1 ) ■ 6q n -x 

T'di^a.qb) - h—{q b ) ) dq b ■ {Sq n - 2 ,Sq n -i) 

Oq J qa=q n -2, 9b=9n-l 

T d (F n ~ 2 v ) - hV'(F n ~ 2 vo)] dq b ■ F n -\5u a - T' d {u Q )dq a ■ 5v 



eU- 


2 »0 


. F n ~ 2 


*5vq + 


K 


■ 5v a 


(f u 


-2* 


8+ ) 

/ v 


8v Q + 


K 


5vq 



for the indicated two-forms 9 + and 6 . Since d(hL d ) = 9 + + 8 , d 2 (hL d ) = = dd + + d6 . Moreover 
the initial conditions vq are arbitrary, hence 



cf S;d = = (F n - Z ) d9 + + d0~ = -(F n - Z ) d9~ + d9 



s< ) 



d$- = (F n - 2 )*d0-. 

Since n is arbitrary, the discrete flow F preserves the symplectic form d0~ . Using backwards error 
analysis, it can be shown that this geometric property guarantees that integrating with F introduces 
no energy drift for a number of steps exponential in h [B] , a highly desirable property when simulating 
molecular dynamic or other Hamiltonian systems whose qualitative behavior is substantially affected 
by errors in energy. 

4. Asynchronous Variational Integrators 



In Section 3.2 an action functional ffify was formulated as the integration of a single discrete La- 



grangian over a single time step size h. Such a construction is cumbersome when modeling multiple 
potentials of varying stiffnesses acting on different parts of the system: to prevent instability one must 
integrate the entire system at the resolution of the stiffest force. Asynchronous variational integrators 
(AVIs), introduced by Lew et al. Q], are a family of numerical integrators, derived from a discrete 
Hamilton's principle, that support integrating potentials at different time steps. Their formulation 
assumes a spatial partition, with each potential depending only on the configuration of a single ele- 
ment; in this exposition, the general arguments set forth by Lew et al. are followed, but the notation 
and derivation departs from their work as necessary to support potentials with arbitrary, possibly 
non-disjoint spatial stencil. 

Let {V 1 } be potentials with time steps h % . Each potential V 1 is concerned with certain moments 
in time — namely, integer multiples of h l — and these moments are inconsistent across triangles. Time 
is therefore subdivided in a way compatible with all triangles: for a r-length interval of time, the set 



S(t) is defined by 

3 M = U U Jh*. 

V* 3=0 

That is, S(t) is the set of all integer multiples less than t of all time steps. S can be ordered, and in 
particular let £(i) be the (i + l)-st least element of S. For ease of notation, also let u> l (j) = C _1 (j^)j 
that is, u) converts the jth timestep of potential i into a global time. 

If n is the cardinality of H, a trajectory of duration t is then discretized by linearly interpolating 
intermediate configurations qo,Qi, ■ ■ ■ ,q n —i, where qi is the configuration of the system at time £(z). 
Velocity is discretized as <7fc+i/2 = t(k+i)-V(k) on ^e se g men t of the trajectory between q^ and qu+i- 
A global action functional of these trajectories is needed, and can be constructed in the natural way: 

n-2 Lt/Zi'J 

SM = E KO' + *) " ^0')] r d fe, Qj+i, £00, £(.? + 1)] - E E ^ '(««'(*)). 

where, for T(g) the kinetic energy of the entire configuration, Td(q a , qb, t a ,tb) = T ( f^Ef^ ) • For use in 

the following, also let T^q a ,q b ,t a ,t b ) = % (» 

No attempt has been made to define a Lagrangian pairing the kinetic and potential energy terms; 
it will be seen that an action defined in this way still leads to a multisymplectic numeric integrator. 
To this end, Hamilton's principle dS g (q) ■ <5q = is imposed for variations <5q = {5qo, . . . , 5q n -i} with 
Sq = 8q n -i = 0. Then S g can be rewritten as 

n— 2 n— 1 

5 fl (q) - E Ktf + !) - to')] r d fo, fc+i, £(j), £(.? + 1)] - E E ^fe), (8) 

where the notation /i. l |£(j) is abused to mean "all indices 2 for which h l evenly divides £(j)," so that 

n-2 n-1 „ 

rf5 9 ( q ) • s q = e ^ [«&, <&+i, e(j). to + 1)] • (<%+i - <%) - E E fti -97 fe) • *& 

= 7d[g„_ 2 ,g„_i,^(n-2),^(n- 1)] -^ n _i - 7^ [go, Qi, C(0), €(!)] • ^0 
V- ui dVi < ^ x 

- 2.^ h 7j~ iin-i) ■ °q n -i 

h'K(n-l) 9 

+ E [ T d fe-n <k> to - i), to)] - id fe. Qj+uZUUU + 1)] - E ^-^-te) ' ^ 

"-2 / ayi \ 

= E \ T 'd fe-i> <?,•, to - *)> to')] - r d te,^+i,€(i),^(i + 1)] - E ^^fe) • <%• 



The Euler-Lagrange equations are then 

dT ,- s dT c s ST ui dV \ ^ ™ 

-gr(qk+i/2) ~ -gr{lk-i/2) = ~ 2^ h ~d^^ k ^ ( 9 ) 

These equations are similar to those derived for synchronous variational integrators fcfh, except that 
only a subset of potentials VJ contribute during each time step. As in the synchronous case, if, as is 
typical, Td(q) is quadratic in q, the system Q gives rise to an explicit numerical integrator that is 
particularly easy to implement in practice. Algorithm [I] gives pseudocode for such integration when 
T ( [ = q T Mq for a mass matrix M; Lew et al. [35] discuss the algorithm in greater detail. 

Algorithm 1 An algorithm for integrating the trajectory given by the AVI Euler-Lagrange equations 
Q adapted from Lew et al. 135 

Let events be (potential, time step, time) triplets E — (V,h,t). 
Denote by qy the configuration subspace on which V depends. 
Let PQ be a priority queue of events, sorted by event times E.t. 
T g <— {T g maintains the value of the simulation clock} 
q <— qo {Set up initial conditions} 

for all Vi do 

Ei <— (Vi, h l , h 1 ) {Add all potentials to the queue as events} 

PQ.push(^) 
end for 
loop 

(V, h, t) «- PQ.pop 

qi-q + (t- T g )q 

4v *— Qv — hMy 1 ^- {Update only those elements affected by this event.} 

_PQ.push(V, h,t + h) {Return the event to the queue, with a new, later time} 

Tg «— t {Update the simulation clock} 
end loop 



4-.1- Multisymplecticity 

The right hand side of ^ depends on £(fc), and so unlike ffih, the Euler-Lagrange equations for AVIs 
are time dependent, and do not give rise to a uniform update rule F(qi-\,qi) i-> [quqi+i)- Instead, 
consider the total, time-dependent flow F k (q a ,qi) n- (qk,qk+i)- Once again, trajectories satisfying 
([9]) are parametrized by vq = (qo,qi), and first variations by 5vo — (8qo,6qi). By restricting to such 
trajectories and variations, the action (pi) can be rewritten as 



[t/VJ-i 

E 

3=0 V J=0 



SiAvi = 'Z[t(j + i)-mT d (p(i>oUtiUti + i)) -E E h i v$(p ,i v +i \ vo )). 



10 



Then, for 



Vf(q a ,q b ) = %(q), 



dS- lA vi{v) ■ &v = dS g (q) ■ 5q 

= T' d [q n - 2 ,q n -i,€(n - 2),£(n - 1)] • 5q n -i - T' d [go, 91,^(0)^(1)] • 6q 



E E *?£(*.-> 



V* h^Un-l) 



dqi 



'in-1 



t; 



^[F n - 2 (v ),an-2),t(n-l) 

J2 E vvi' \F n -\v a ) 



■5q n _ x -T' d [^,^(0)^(1)] -Sq 
Sq n -i 



I'O 



•' - 2 *6 + ) V0 ■ 5v Q 
for one-forms 6~ and 9 + . Once again 



= e~ n ■ 6v + ff^_ ■ F" 



F r ' 



= d 2 5 iA vi = dJS~ + F n - 2 *de + , 



(10) 



but unlike when the action was a sum of Lagrangians, from the multisymplectic form formula ( 10 1 



there is no way of relating d0~ to d9 + , and thus discrete symplectic structure preservation is not 
recovered. Nevertheless, Lew et al. [1] conjecture that this multisymplectic structure leads to the good 
energy behavior observed for AVIs. 



5. Discrete Penalty Layers 

The above reformulation of AVIs can be leveraged to resolve collisions with guaranteed perfect 
robustness, and via momentum-symplectic integration, so that the energy behavior of the system as a 
whole remains good. Consider a standard penalty force approach, which for every two elements A, B 
and surface thickness 77 defines the gap function 

g n (q) = inf ||a-6|| - 2rj 

a^A.b^B 

measuring the proximity of A to B. 

The penalty potential is then defined as 



V{q) 



g v {q) > 

kg v (q) 2 g{q) < 0, 



where k is a user-specified stiffness. As previously discussed, V alone does not robustly prevent 
interpenetrations: the potential can be viewed as placing a spring between the approaching elements, 
and for sufficiently large relative momentum in the normal direction, the spring will fully compress, 
then fail. However, consider placing an infinite family of potentials Vi, I = 1,2,..., between the 
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Figure 1: Plots of the potential energy of the first three layers as a function of gap function g (left), and a plot of 
the total potential energy contributed by all layers < n for n = 1,2,3 (right). Notice the potential energy diverges as 
separation distance approaches 0, guaranteeing that collision response is robust. 



primitives, where 



Vi(q) 



l 3 kg 2 



n/i 



9r,/l{q) > o 

g v /M) < o. 



The region ^_. < d(q) < ^, where exactly n of the potentials are active, is called the n-th discrete 
penalty layer. Figure [T] shows a plot of the potential energy of the first few potentials for the case 
T] = k = 1, as well as the cumulative potential energy of all of the potentials. 
The total potential energy of the springs when fully compressed is 



OO „ oo 



i=i 

which diverges. The infinite array of potentials is guaranteed to stop all collisions. This guarantee in 
no way depends on the chosen stiffness k: although performance and error will vary with the choice 
of stiffness, unlike for penalty forces the stiffness does not affect the guarantee. The method is always 
guaranteed to be robust. 

There is one obstruction to implementing this scheme in practice: integrating the Z-th spring stably 
and with good energy behavior requires a time step proportional to j^, which vanishes as / — > oo. 
Using a traditional integrator, one could decide ahead of time to only simulate the first few springs — 
but then the guarantee that no penetrations will occur is lost, and the simulation must be run at a 
prohibitively small time step. AVfs, with the above modifications, and a bit of extra bookkeeping, are 
a first step towards alleviating the problem, by allowing the user to assign each spring its own time 
step. This bookkeeping is now described, in terms of modifications to the basic Algorithm [T] 



6. The Asynchronous Algorithm 

AVIs allow each penalty layer to be assigned a different time step, so that less stiff (/ small) 
layers can take large time steps regardless of the presence of the stiffer layers. However, it is still not 
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possible as a practical matter to integrate the system, since arbitrarily large I would need arbitrarily 
small time steps, and the global time in Algorithm [l] would never advance. The following observation 
surmounts this obstacle: at any time during a well-posed simulation, the number of layers that are 
exerting a non-zero force, or that are active, is finite. More precisely, a simulation is well-posed if 
its total energy over time is bounded — that is, if the simulation begins in a non-penetrating state; all 
prescribed, infinite-mass bodies are stationary; and only a finite amount of energy is added over time 
in the form of external forcing. Inactive penalty potentials can be ignored by Algorithm [l] entirely, 
since they do not change configurational velocity, and the position integration that would take place 
during the handling of an inactive potential can just as well be done by the following event. Therefore 
the simulation would be guaranteed to never stop making progress if there is a lower bound for the 
amount of global time T g that elapses with the processing of any event. Such a lower bound exists 
if there is a way to detect which penalty potentials are active or inactive at all times and remove all 
inactive events from the priority queue PQ. 

Suppose that at the start of the simulation, all penalty layers are inactive. Thus no penalty layer 
events are needed on the queue. For each pair of simulation elements, the time t a that the first penalty 
layer would become active (assuming all elements continue along the trajectory described by their 
initial velocities) can be calculated, and the corresponding event added to the queue at that time. 
Such an approach suffers from two problems, however. Firstly, solving for the time when the gap 
function will be zero is easy in some cases, such as if the elements are two spheres or two planes, 
but can involve expensive root solves in others, such as if the elements are two non-rigid triangular 
elements of a thin shell simulation. Secondly, the times computed are fragile: should any event alter 
the velocity of one of the elements (such as a material force, or gravity, or another penalty force if 
one of the elements collides with a third party) the activation time is no longer valid and must be 
recomputed. 

Instead of an exact time, only a conservative guarantee, or certificate |17j . that the first penalty 
layer will not be active before some time t c (where necessarily t c < t a ) is truly needed. For example, 
one certificate is the existence of an 2?7-thick planar slab S that separates the two elements up until time 
t c , where 77 is the thickness of the first penalty layer. For an m-dimensional configuration space, such 
a planar slab is understood to be an extrusion of an (to — 1) — dimensional affine subspace. Concretely, 
let w be a unit vector in W 11 , im be m — 1 linearly independent vectors in K m orthogonal to w, and p 
a point in K m . Then the slab S wp is the set 



S w , p = \p + aw + J^ fiiWi - r\ < a < 77, j3 % € R j 



If such a slab separates the two elements, the first penalty layer cannot become active before t c . This 
certificate can be placed as an event on the queue, with time t c . The certificate might then suffer 
several fates: [T5] 

• An event modifies the velocity of one of the elements before time t c . The certificate placed on 
the queue is then no longer valid until time t c , but instead until a new time t' c which may be 
sooner or later than t c . The algorithm must thus reschedule the certificate, by removing its event 
from the queue, and reinserting it at the appropriate new time. 

• The certificate event is popped from the queue without incident, but it is possible and convenient 
to find a new separating slab that guarantees the penalty layer does not activate before time 
t' c > t c . This new certificate can then be pushed on the queue for time t' c . 
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• The certificate event is popped from the queue without incident, but finding a new slab is 
impossible, costly, or a slab can be found, but the new time t' c is judged heuristically to be too 
near t c . The first penalty layer may then be activated early: doing so affects the efficiency, but 
not the correctness, of the simulation. Simultaneously, the algorithm searches for an ?7-thick 
separating slab to serve as a certificate that layer two is not yet active, and the whole process 
described above is repeated. 

Detecting when a penalty layer event becomes inactive, and should be removed from the queue, 
is much simpler than detecting layer activation: whenever a penalty force for layer n is integrated, 
the algorithm simply checks if the force applied was 0. If so, and if the two elements in question are 
separating, layer n is now inactive: it is not pushed back onto the queue (and instead a separating 
slab of thickness -q/n is sought.) 

It is very important to note that when an event becomes active and is added back into the event 
priority queue, it is done so at a time that is an integer multiple of its timestep from its last time 
of integration. That is, those times when integration would do nothing have been optimized away, 
but the potential's "integration clock" has not been tampered with or realigned, since every potential 
having a fixed-size time step was fundamental to the proof that asynchronous variational integration 
is multisymplectic. The spring-on-a-plane example described below underlines the danger of failing to 
maintain such a fixed time step. 

For an event E, denote all simulation elements on which E depends the support of V. Denote 
all simulation elements whose velocities are modified by E the stencil of E. For force integration 
events, there is no distinction between stencil and support. Certificates have a support, but no stencil. 
Algorithm [2] uses this terminology to incorporate the above into the AVI algorithm. 

Algorithm 2 Proposed algorithm for asynchronous contact resolution. 
Let force events be (potential, time step, time) triplets E — (V, h,t). 
Let PQ be a priority queue of events, sorted by event times E.t. 
T g <— {T g maintains the value of the simulation clock} 
q <— <7o {Set up initial conditions} 

q <- go 

Push non-penalty (e.g. material) events on the queue 

for all pairs of elements K\ , K2 do 

E <- FindCertificate^i, K 2 ) 

PQ.push{E) 
end for 
loop 

E <- PQ.pop 

q <- q + (E.t - T g )q 

if E is a force event then 
handleForceEvent(PQ, E) 

else 

handlcCcrtificateEvent(PQ, E) 

end if 

T g 4— E.t {Update the simulation clock} 
end loop 
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Algorithm 3 handlcForceEvent 



Require: Priority queue of events PQ and force event E that needs processing 

{Processing a force event E is a three-step process: integrating the force, rescheduling all events 
whose support depends on E's stencil, and lastly, resccduling E itself.} 
for all i in Stencil(-E) do 

Q% <— Q% — (E.h)M^ l d Q q v {Update only those elements affected by this event.} 
end for 

{Reschedule all events whose support depends on E's stencil} 
for all Certificate events E' with Stencil(-E) n Support(S') ^ do 

P<5.remove(-E') 

Schedule^') 

PQ.push(£') 
end for 

{If E was a penalty force event, it exerted force, and the two primitives in question are separating, 
then we no longer need it} 
if E is a penalty force event and 3 P = then 

if E.V.K1 and E.V.K2 have positive relative velocity (are separating) then 
return 

end if 

if addCertfficate(.E.V.Jn,.E.V..K'2) then 
return 

end if 
end if 

{Otherwise, reschedule E itself} 
PQ.push(V, h,t + h) 



Algorithm 4 handlcCertificatcEvcnt 



Require: Priority queue of events PQ and certificate event E that needs rescheduling 
if not addCertificate^.in, E.K2) then 

{Finding a new certificate failed. We must thus activate a penalty force, one layer deeper than 
the deepest currently active penalty force event.} 
CurLayer <— max{ pcnalty cvonts E > on qucuo f or e.k\ and e.K2} E -layer 
E' <— new PenaltyForceEvent(i?.i ; £'l, E.K2, CurLayer + 1) 
-PQ.push(_E') {Push the appropriate penalty force event on the queue} 
end if 
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Algorithm 5 addCertificate 

Require: Priority queue of events PQ, and two elements K\ and K2 

{Attempts to find a certificate for the collision of K\ against K1 and add it to the queue. Returns 

true if one was found.} 

E' <- FindCertificate^l, K2) 

if E' was successfully found then 
PQ.push(FindCertificate(!n, K2)) 
return true 

end if 

return false 



In Algorithm [2] and its accompanying subalgorithms, the behavior of the functions FindCertificate 
and Schedule will depend on the type of certificate chosen. FindCertificate returns a new certificate 
for a given pair of elements, if possible and practical, and Schedule computes the time a certificate 
becomes invalid, as described in the paragraphs above. For thin shell simulation, where all simulation 
elements are convex triangles, edges, and vertices, separating slabs serve as ideal certificates, since it 
is cheap to compute Schedule, in this case by calculating element-plane intersection times. Although 
any choice of certificate, and heuristic for when to abort searching for a new certificate, preserves the 
correctness of the algorithm, the progress property described in the first paragraph of this section relies 
on the certificates efficiently weeding out inactive events so that some certificate is found before all 
(infinitely many) layers for a pair of elements are activated. No problems have been observed using 
separating slabs for thin-shell simulations, but different certificates may be needed, e.g., for concave 
rigid bodies. 

6.1. Further Optimizations 

The technique explored in the previous section, of finding a sequence of conservative certificates 
guaranteeing that some property holds, instead of calculating an exact time when that property stops 
holding, is the central idea behind a wide class of algorithms known as Kinetic Data Structures (KDSs) 
|17j . In the case described above, the property was inactivity of a given penalty layer. KDSs are 
particularly well-suited for an asynchronous approach, since certificate expiration times may not all 
align to some convenient simulation clock, and the required rescheduling of certificates/searching for 
new certificates can reuse the priority queue data structure already needed for force integration events. 
To improve the efficiency of the implementation used to create the examples below, several more KDSs 
in addition to the separating slabs discussed above were implemented: a bounding volume hierarchy 
[39] was used to take advantage of the fact that spatially distant elements are unlikely to collide, 
separation lists [?D] to optimize the bookkeeping of this hierarchy, and a novel KDS was devised to 
leverage the observation that high-frequency, low amplitude oscillations in velocity do not significantly 
change a separation slab's expiration time, so that rescheduling is in many cases unnecessary. All of 
the improvements are described in greater detail in |16j . 

7. Dissipation 

The framework, as described so far, gives near-perfect long-time energy conservation. In the real 
world, however, many dissipative phenomena are observed — for instance friction, spring damping, 
and non-unit coefficients of restitution during collisions. Several simple modifications can be made 
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to the proposed method to take such dissipation into account. Qualitatively, these have performed 
well in practice: energy seems to behave well over long times for dissipative systems analogously to 
the near- conservation observed for Hamiltonian systems, but a theoretical understanding of this good 
behavior remains future work. 

7.1. Coefficient of Restitution 

It is often desirable to simulate semi-elastic or inelastic collisions. A simple modification to the 
potential Vi allows the use of arbitrary coefficients of restitution e: 



V l { q) = [ Q thrt9)><> 

\l 3 ksg 7l/ i(q) 2 g v /i(q)<0, 



where s is e if the primitives are separating, 1 otherwise. The penalty layers exert their full force 
during compression, then weaken according to the coefficient of restitution during decompression. 

This approach, while simple, does have a limitation in the inelastic limit e = 0: due to error 
introduced by numerical integration, two colliding primitives may have non-zero, though small, post- 
response normal relative velocity. The magnitude of this velocity is at most kr//l, so it can be limited 
by choosing a small enough base stiffness k. 

7.2. Friction 

The Coulomb friction model is a simple approximation to kinetic friction: at a point of contact 
between two bodies, the Coulomb force has magnitude /u|-F n |, where /x is a coefficient of friction and 
F n is the normal force at the contact points, and has direction opposite the relative tangential motion 
of the contact points. 

Whenever an impulse is applied during integration of a penalty layer, a corresponding frictional im- 
pulse can also be applied. Just as increasingly stiff penalty forces are applied for contact forces, friction 
forces are increasingly applied (equal to /u|-F n |) to correctly halt high-speed tangential motion. Notice 
that these friction forces, like the material and contact penalty forces, are applied asynchronously: 
every layer applies friction independently at its own time step. 

This simple, asynchronous formulation of friction fits very naturally into the framework of AVIs. 
Unfortunately, it is unsuitable for simulations featuring static friction, such as a block of wood resting 
on an inclined plane. The above formulation, with friction applied piecemeal during penalty inte- 
gration, is reactive instead of proactive, and in simulations of this type the block of wood has been 
observed "creeping" down the incline no matter how high a coefficient of restitution is chosen. A more 
comprehensive model of friction compatible with the AVI framework, which correctly handles static 
friction, remains future work. 

8. Results 

8.1. Spring on a plane 

As a simplest didactic demonstration of the proposed method, three experiments were conducted. 
A vertical spring of unit rest length, stiffness, and endpoint masses began each of the three simulations 
stationary a unit height above a fixed horizontal plane. The springs fell under a gravitational force of 
strength lm/s 2 , with impact handled in one of three different ways: 

In the first experiment, gravity and the stretching force were integrated synchronously, and an 
instantaneous impulse was applied whenever the bottom of the spring touched the plane. Figure [2] 
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Figure 2: The total energy of a spring bouncing on a plane integrated by mixing variational material force integration 
with impulse-based collision response (blue), the energy of this simulation when the alignment of the integration clocks is 
not respected (maroon), and the energy of this simulation when using the proposed method (brown). The exact energy 
is shown as a dashed black line. 

shows energy over time when using this method, compared to expected perfect energy conservation. 
Energy in this experiment, far from being conserved, took a random walk. 

In the second experiment, all forces were integrated asynchronously using the proposed method. 
The thickness r\ was chosen to be 0.1, and the penalty base stiffness k, 1000. Energy in this case was 
well-conserved over long time: although energy experiences high-frequency, low-amplitude oscillations, 
there was no drift. 

The importance of respecting the integrity of each potential's integration clock is highlighted in 
the third experiment. Instead of adding a force event onto the priority queue at an integer multiple of 
its time step, the events are added at the precise moment when each layer becomes active. As can be 
seen from the resulting energy plot (Figure K|), energy is no longer well-conserved, but instead seems 
to increase monotonically over time. 

8.2. Balls of Particles 
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Figure 3: A rigid box containing 900 spheres with random initial velocity, several minutes after the start of the simulation. 



18 



200 r 



ft 150 



E 100 

o> 
E 
o 
2 50 



6000 r 
5000 
4000 
i 3000 
2000 
1000 



10 



20 30 

Time ( s ) 



40 



50 



10 15 

Time ( s ) 



20 



25 



Figure 4: Left: The energy (brown) and momentum (orange) of the box of the spheres as a function of time. Right: The 
energy of the box, with gravity added, over time. 

As an example that involves more collisions, consider a fixed 3 m x 3 m square box. Inside this box 
900 spheres of radius 10 cm were uniformly distributed, each of which was given a random velocity of 
magnitude between and lOm/s. Figure [3] depicts this box after several minutes have elapsed. Energy 
over time is plotted in Figure El (left), and it is again almost perfectly conserved. Figure 13] (center) 
plots the magnitude of total momentum of the box over time, and it is exactly conserved, as expected 
since a multisymplcctic integrator is used. Gravity (9.8 m/s 2 ) was added to the box and again energy 
plotted over time (Figure |4J right), and again good energy behavior was observed. 




Figure 5: The energy of the box of 900 spheres under different coefficients of restitution: from top to bottom, 1.0, 0.9, 
0.8, 0.7, 0.5, 0.2, 0.0. 

As a test of controllable dissipation by using a coefficient of restitution, the box with gravity was 
resimulatcd several times using different coefficients of restitution. Figure [5] shows the resulting energy 
plots. For any chosen coefficient of restitution, the non-conservative energy behavior is qualitatively 
as one would expect. 

8.3. Sphere-Plate Impact 

The impact experiment of a spherical shell against a thin plate, as described in Cirak and West's 
article on Decomposition Contact Response (DCR) [33] . was reproduced using the proposed framework. 
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A sphere of radius 12.5 cm approaches a plate of radius 35 cm with relative velocity 100 m/s. Both 
the sphere and the plate have thickness 0.35 cm. The time steps of the material forces (stretching and 
bending) are 10 _7 s (the same as those chosen by Cirak and West.) 
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Figure 6: Total energy over time of a thin sphere colliding against a thin plate, simulated using the proposed contact 
response method (right) compared to data provided for decomposition contact response I33| (left). 

Figure [6] compares energy over time when this simulation is run using both the proposed method 
and DCR. Using the former there is no noticeable long-term drift; closely examining the energy data 
reveals the high-frequency, low-amplitude, qualitatively-negligible oscillations characteristic of sym- 
plectic integration. The latter introduces noticeable artifical energy decay. 

8.4- Large-scale Three-dimensional Examples 

Harmon et al. [16] describe a series of optimizations that improve the efficiency of Algorithm [2j 
These optimizations were incorporated to form our Asynchronous Contact Mechanics (ACM) code. 
This code continues to yield qualitatively good results when scaled to additional large-scale problems. 




Figure 7: Simulated tying of a cloth reef knot (left) and bowline knot (right). 

Two thin rectangular 27 cm x 2 cm ribbons were modeled as thin shells of 5321 vertices, subject to 
constant-strain triangle stretching forces |JT] (stiffness 750) and discrete shell bending forces formulated 
by Grinspun et al. [H] (stiffness 0.05). These ribbons were positioned into a loose reef knot by an 
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artist. The knot was then tightened by constraining the end of the ribbon to move apart at lOcm/s, 
and running the simulation. 

Figure[7J left, shows the ribbon after 2 seconds. Since the velocities of the ends of the ribbons were 
constrained, the knot material became arbitrarily stretched once the knot was tight. The forces pressing 
the two ribbons into each other thus grew unbounded, but the two ribbons never interpenetrated, nor 
were other collision-related artifacts observed. It should be stressed that this good behavior did not 
require the tweaking of the penalty stiffnesses nor any other artificial parameters. 

As a second large-scale example, a ribbon similar to the ones in the reef knot simulation was 
positioned by an artist into a loose bowline knot tied around a cylindrical thin shell of 1334 vertices. 
The bowline was then tightened by fixing one end of the ribbon and constraining the other to move 
away from the cylinder at 10 cm/s. Again, the knot successfully tightened with no penetrations or 
other artifacts (figure [7j right). 

8.5. Sphere and Wedge 





Figure 8: A sphere falling into a wedge, at the beginning of the simulation (left and center) and 0.42 seconds later, after 
the sphere has reflected off of the wedge (right). The center figure shows the mesh elements of the bodies. 

Inspired by Pandolfi et al. [lOj . a rigid thin-shell sphere was dropped into a wedged formed by two 
thin shell triangular prism, shown in Figure [8J Each prism has an isosceles base with width 12.92 cm 
and height 20.05 cm, and length 38.41 cm. The prisms contain 71 vertices each. The sphere contains 
92 vertices, has radius 4.97 cm and begins the simulation 20.84 cm above the ground plane on which 
the prisms rest. The sphere has initial downwards velocity of —100 cm/s (no gravity). The sphere and 
shells use the same thin shell model as the debris in the above trash compactor example, with bending 
and stretching stiffness parameters 100000 and 50000 respectively. 

As the sphere descends, it enters into multiple contact with the faces of the wedge, which undergo 
elastic deformation and high-frequency vibration. Despite the large areas of simultaneous contact and 
high velocity at the time of impact, the energy of this system, plotted in Figure |9j exhibits good 
behavior and does not drift. 



8.6. Draping on Spikes 

ACM's ability to robustly handle degenerate geometry was tested by dropping two 1994-vertex, 
15 cm x 50 cm cloth meshes (stretching stiffnes 500, bending stiffness 0.1, stretching damping 1.0, 
bending damping 0.1) on top of a rigid 20 cm x 20 cm plate from which protrude 36 7.8 cm spikes 
(see figure 10). Each spike was modeled using six highly-degenerate, sliver triangles: each triangle's 
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Figure 9: The energy of the wedgc-sphcrc system as a function of time. The energy docs not drift over the course of the 
simulation. 




Figure 10: Two cloth rectangles were draped on a bed of spikes. The system at the start of the simulation (left), and 
after the cloth has come to rest (right). 




Figure 11: A close-up of one of the spikes; the spike has been rotated clockwise 90 degrees to conserve space. Each spike 
is composed of six triangles with apex angle 3.47 degrees. 



most acute angle measures 3.47 degrees (see figure 11). The cloth was allowed to fall under gravity 
(9.8 m/s 2 ) and drape on top of the spikes until it had come to rest, 
other artifacts were observed. 



No penetrations, oscillations, or 

After the cloth came to rest, the bottom cloth was pulled out from under the top one by constraining 

The 
No 



one side of the cloth to move at 10 cm/s parallel to and away from spiked plate; see figure 12 



bottom cloth scraped against the spikes and slid, with no dissipation, against the top cloth. 
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Figure 12: After the cloth came to rest, the bottom cloth was pulled out from under the top one. The simulation 3 
seconds after pulling began. 



interpenetrations occured. 
8.7. Trash Compactor 





Figure 13: Walls close in and compress various thin-shell objects. The beginning (left) and end (right) of the simulation. 

Various coarse thin-shell solid objects (platonic solids, tori, etc.) modeled as triangle meshes were 
placed in a rectangular box measuring 71.5 cm x 36.7 cm x 9.3 cm. The four sides were scripted to close 
in and compress the objects within: the length at 20 cm/s, and the width at 10 cm/s. All objects were 
given the same material parameters (stretching stiffness 1000, bending stiffness 10, stretching damping 
15, bending damping 0.5) and held to the floor of the box by gravity (9.8 m/s 2 ). Figure 
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shows the 
box at the beginning of the simulation, and after the simulation had run for 3.4 seconds. A simple 
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plastic deformation model, described by Bergou et al. [33], allowed the objects to crush plastically 
when stressed by the encroaching walls. Nevertheless, the material forces acting on the objects grew 
larger as the box decreases to a small fraction of its original volume, yet no object penetrated any 
other object or wall, as guaranteed by the method. 

9. Conclusion and Future Work 

A framework for asynchronous, structure-preserving handling of contact and impact has been pre- 
sented. Provable guarantees were established for this framework: impact handling is robust, allowing no 
penetrations or tunneling; the good properties of AVIs are preserved, such as a discrete Noether's The- 
orem and discrete multisymplectic structure; and for well-posed problems, the amount of computation 
required to simulate the problem is bounded and in particular finite. Good long-time energy behavior, 
conjectured to accompany multisymplectic structure, was confirmed empirically by both didactic and 
challenging, large-scale experiments. Modifications to allow for simple dissipative phenomena, such as 
a coefficient of restitution, were described. Data structures and algorithms to improve performance, 
such as the use of separating slabs to prune inactive penalty layers, were briefly discussed. Imple- 
mentation details for these and other optimizations, as well as ideas for future improvements to the 
algorithms that promise to substantially decrease computation time, are treated more comprehensively 
by Harmon et al. |16j . 

Missing from the basic asynchronous contact framework described by this paper is comprehensive 
handling of friction, particularly static friction. Static friction conflicts fundamentally with asynchrony: 
in an asynchronous simulation, contact between a pair of elements is resolved piecemeal, by summing 
the impulses at many different times contributed by many different penalty layers. At any given 
moment of time it is unclear how to define a total normal force, an element necessary for the robust 
treatment of even the most elementary static friction models. Successfully merging the handling of 
friction with the asynchronous framework, to allow simulations of systems such as a standing house of 
cards, remains a challenging area for future research. 
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